"""
GCR.jl Integration Layer

Provides seamless integration with the GCR.jl framework for solving Poisson equations
in the lid-driven cavity solver. Handles both pure Julia and PETSc backends with
appropriate preconditioners and solver options.
"""

# Helper functions for GCR result compatibility
"""
    get_solution_vector(result)

Extract solution vector from GCR result, handling different result types.
"""
function get_solution_vector(result)
    if hasfield(typeof(result), :x)
        return result.x
    elseif isa(result, AbstractVector)
        return result
    else
        error("Unable to extract solution vector from result of type $(typeof(result))")
    end
end

"""
    get_convergence_status(result)

Extract convergence status from GCR result.
"""
function get_convergence_status(result)
    if hasfield(typeof(result), :converged)
        return result.converged
    else
        # Assume converged if we got a result
        return true
    end
end

"""
    get_iteration_count(result)

Extract iteration count from GCR result.
"""
function get_iteration_count(result)
    if hasfield(typeof(result), :iterations)
        return result.iterations
    else
        return 0  # Unknown iteration count
    end
end

"""
    get_residual_norm(result)

Extract final residual norm from GCR result.
"""
function get_residual_norm(result)
    if hasfield(typeof(result), :residual_norm)
        return result.residual_norm
    else
        return NaN  # Unknown residual norm
    end
end

"""
    PoissonSolver

Abstract type for Poisson equation solvers.
"""
abstract type PoissonSolver end

"""
    GCRPoissonSolver <: PoissonSolver

Poisson solver using GCR.jl framework with configurable backend.

# Fields
- `matrix::SparseMatrixCSC{Float64, Int}`: Assembled Poisson matrix
- `id_map::Array{Int}`: Index mapping from full grid to interior nodes
- `interior_indices::Vector{CartesianIndex}`: Indices of interior nodes
- `solver_options::NamedTuple`: GCR solver configuration
- `backend::Symbol`: Backend type (:julia, :petsc, :gpu)
- `preconditioner::Any`: Preconditioner object
"""
mutable struct GCRPoissonSolver <: PoissonSolver
    matrix::SparseMatrixCSC{Float64, Int}
    id_map::Array{Int}
    interior_indices::Vector{CartesianIndex}
    solver_options::NamedTuple
    backend::Symbol
    preconditioner::Any
    
    function GCRPoissonSolver(matrix, id_map, interior_indices; 
                             backend=:julia, 
                             tol=1e-8,
                             maxiter=1000,
                             restart=30,
                             preconditioner=:diagonal,
                             verbose=false)
        
        # Create solver options
        options = (
            tol = tol,
            maxiter = maxiter,
            restart = restart,
            verbose = verbose,
            save_history = true
        )
        
        # Initialize preconditioner
        P = setup_preconditioner(matrix, preconditioner, backend)
        
        return new(matrix, id_map, interior_indices, options, backend, P)
    end
end

"""
    AdaptivePoissonSolver <: PoissonSolver

Adaptive Poisson solver that automatically selects the best preconditioner
and solver configuration based on problem properties and runtime performance.

# Fields
- `matrix::SparseMatrixCSC{Float64, Int}`: Assembled Poisson matrix
- `id_map::Array{Int}`: Index mapping
- `interior_indices::Vector{CartesianIndex}`: Interior node indices
- `solver_options::NamedTuple`: Base solver configuration
- `backend::Symbol`: Backend type
- `preconditioner_candidates::Vector{Symbol}`: Available preconditioners
- `best_preconditioner::Symbol`: Currently best performing preconditioner
- `performance_history::Dict`: Performance tracking
- `refinement_enabled::Bool`: Whether iterative refinement is enabled
"""
mutable struct AdaptivePoissonSolver <: PoissonSolver
    matrix::SparseMatrixCSC{Float64, Int}
    id_map::Array{Int}
    interior_indices::Vector{CartesianIndex}
    solver_options::NamedTuple
    backend::Symbol
    preconditioner_candidates::Vector{Symbol}
    best_preconditioner::Symbol
    performance_history::Dict{Symbol, Vector{Float64}}
    refinement_enabled::Bool
    current_solver::Union{PoissonSolver, Nothing}
    
    function AdaptivePoissonSolver(matrix, id_map, interior_indices;
                                  backend=:julia,
                                  tol=1e-8,
                                  maxiter=1000,
                                  restart=30,
                                  verbose=false,
                                  refinement_enabled=true)
        
        # Determine available preconditioners based on problem size and backend
        candidates = determine_preconditioner_candidates(matrix, backend)
        
        # Create base options
        options = (
            tol = tol,
            maxiter = maxiter,
            restart = restart,
            verbose = verbose,
            save_history = true,
            adaptive_tolerance = true
        )
        
        # Initialize performance tracking
        history = Dict{Symbol, Vector{Float64}}()
        for candidate in candidates
            history[candidate] = Float64[]
        end
        
        new(matrix, id_map, interior_indices, options, backend,
            candidates, candidates[1], history, refinement_enabled, nothing)
    end
end

"""
    FallbackPoissonSolver <: PoissonSolver

Fallback solver using Julia's built-in linear algebra when GCR.jl is not available.
"""
mutable struct FallbackPoissonSolver <: PoissonSolver
    matrix::SparseMatrixCSC{Float64, Int}
    id_map::Array{Int}
    interior_indices::Vector{CartesianIndex}
    factorization::Any
    
    function FallbackPoissonSolver(matrix, id_map, interior_indices)
        # Use built-in sparse direct solver or iterative methods
        try
            # Try sparse LU factorization for moderate-sized problems
            F = lu(matrix)
            return new(matrix, id_map, interior_indices, F)
        catch
            # Fall back to unfactorized form (use CG at solve time)
            return new(matrix, id_map, interior_indices, nothing)
        end
    end
end

"""
    create_poisson_solver(domain::LShapeDomain{D}, h::Float64, options::LidDrivenOptions) 
    -> PoissonSolver where D

Create an appropriate Poisson solver based on available packages and options.

# Arguments
- `domain::LShapeDomain{D}`: Domain geometry
- `h::Float64`: Grid spacing
- `options::LidDrivenOptions`: Solver configuration

# Returns
- `solver::PoissonSolver`: Configured Poisson solver
"""
function create_poisson_solver(domain::LShapeDomain{D}, h::Float64, options::LidDrivenOptions) where D
    # Assemble Poisson matrix
    A, id_map, interior_indices = assemble_poisson_matrix(domain, h)
    
    if options.verbose
        print_matrix_info(A, domain)
    end
    
    # Choose solver based on availability and options
    if HAS_GCR
        return GCRPoissonSolver(A, id_map, interior_indices,
                               backend = options.solver,
                               tol = options.tol * 1e-3,  # Tighter tolerance for Poisson
                               maxiter = min(1000, 2 * size(A, 1)),
                               preconditioner = options.preconditioner,
                               verbose = options.verbose)
    else
        if options.verbose
            @info "Using fallback Julia solver (GCR.jl not available)"
        end
        return FallbackPoissonSolver(A, id_map, interior_indices)
    end
end

"""
    solve_poisson!(solver::GCRPoissonSolver, rhs_field::Array{Float64}, 
                  solution_field::Array{Float64}) -> NamedTuple

Solve the Poisson equation ∇²ψ = rhs using GCR.jl framework.

# Arguments
- `solver::GCRPoissonSolver`: Configured GCR solver
- `rhs_field::Array{Float64}`: Right-hand side field (full grid)
- `solution_field::Array{Float64}`: Solution field (modified in-place)

# Returns
Named tuple with solve diagnostics:
- `converged::Bool`: Whether solver converged
- `iterations::Int`: Number of iterations
- `residual_norm::Float64`: Final residual norm
- `solve_time::Float64`: Solver time in seconds
"""
function solve_poisson!(solver::GCRPoissonSolver, rhs_field::Array{Float64}, 
                       solution_field::Array{Float64})
    
    # Extract RHS values at interior nodes
    rhs_interior = extract_interior_values(rhs_field, solver.interior_indices)
    
    # Extract initial guess from solution field
    x0 = extract_interior_values(solution_field, solver.interior_indices)
    
    # The matrix A represents -∇², so we need to solve: -A*ψ = -rhs
    # This is equivalent to: A*ψ = rhs
    b = rhs_interior  # Since our matrix assembly already has the correct sign
    
    # Solve using GCR.jl with improved error handling
    t_start = time()
    
    try
        if HAS_GCR && GCR_MODULE !== nothing
            # Use dynamic dispatch to call GCR solver
            solve_gcr_func = getfield(GCR_MODULE, :solve_gcr)
            
            # Prepare solver arguments with proper types
            solver_args = Dict{Symbol, Any}(
                :tol => solver.solver_options.tol,
                :maxiter => solver.solver_options.maxiter,
                :restart => solver.solver_options.restart,
                :verbose => solver.solver_options.verbose,
                :save_history => solver.solver_options.save_history,
                :x0 => x0
            )
            
            # Add preconditioner if available
            if solver.preconditioner !== nothing
                solver_args[:preconditioner] = solver.preconditioner
            end
            
            # Call GCR solver with dynamic dispatch
            result = solve_gcr_func(solver.matrix, b; solver_args...)
            
            # Extract solution and check result structure
            solution_vector = get_solution_vector(result)
            
            # Scatter solution back to full field
            scatter_interior_values!(solution_field, solution_vector, solver.interior_indices)
            
            solve_time = time() - t_start
            
            return (
                converged = get_convergence_status(result),
                iterations = get_iteration_count(result),
                residual_norm = get_residual_norm(result),
                solve_time = solve_time,
                solver_info = "GCR.jl $(solver.backend)",
                gcr_result = result  # Store full result for debugging
            )
        else
            error("GCR.jl not available but GCRPoissonSolver was created")
        end
        
    catch e
        @error "GCR solve failed: $e"
        @debug "Matrix properties: size=$(size(solver.matrix)), nnz=$(nnz(solver.matrix))"
        @debug "RHS properties: length=$(length(b)), norm=$(norm(b))"
        
        # Fall back to emergency solver
        return solve_poisson_fallback!(solver, rhs_interior, solution_field)
    end
end

"""
    solve_poisson!(solver::FallbackPoissonSolver, rhs_field::Array{Float64}, 
                  solution_field::Array{Float64}) -> NamedTuple

Solve Poisson equation using Julia's built-in linear algebra.
"""
function solve_poisson!(solver::FallbackPoissonSolver, rhs_field::Array{Float64}, 
                       solution_field::Array{Float64})
    
    rhs_interior = extract_interior_values(rhs_field, solver.interior_indices)
    b = rhs_interior
    
    t_start = time()
    
    try
        if solver.factorization !== nothing
            # Use pre-computed factorization
            x = solver.factorization \ b
            converged = true
            iterations = 1  # Direct method
            residual_norm = norm(solver.matrix * x - b)
        else
            # Use conjugate gradient
            x0 = extract_interior_values(solution_field, solver.interior_indices)
            x, ch = cg(solver.matrix, b, x0, tol=1e-8, maxiter=min(1000, length(b)))
            
            converged = ch.isconverged
            iterations = ch.iters
            residual_norm = ch.residuals[end]
        end
        
        solve_time = time() - t_start
        
        # Scatter solution back
        scatter_interior_values!(solution_field, x, solver.interior_indices)
        
        return (
            converged = converged,
            iterations = iterations,
            residual_norm = residual_norm,
            solve_time = solve_time,
            solver_info = "Julia fallback"
        )
        
    catch e
        @error "Fallback Poisson solve failed: $e"
        solve_time = time() - t_start
        
        return (
            converged = false,
            iterations = 0,
            residual_norm = Inf,
            solve_time = solve_time,
            solver_info = "Failed"
        )
    end
end

"""
    setup_preconditioner(matrix::SparseMatrixCSC{Float64, Int}, 
                        precond_type::Symbol, backend::Symbol) -> Any

Set up preconditioner for the Poisson matrix.

# Arguments  
- `matrix::SparseMatrixCSC{Float64, Int}`: Poisson matrix
- `precond_type::Symbol`: Preconditioner type (:none, :diagonal, :ilu, :amg)
- `backend::Symbol`: Solver backend (:julia, :petsc, :gpu)

# Returns
- Preconditioner object appropriate for the backend
"""
function setup_preconditioner(matrix::SparseMatrixCSC{Float64, Int}, 
                             precond_type::Symbol, backend::Symbol)
    
    if precond_type == :none
        if HAS_GCR && GCR_MODULE !== nothing
            # Use GCR.jl identity preconditioner
            try
                IdentityPreconditioner = getfield(GCR_MODULE, :IdentityPreconditioner)
                return IdentityPreconditioner()
            catch
                return nothing
            end
        else
            return nothing
        end
        
    elseif precond_type == :diagonal
        if HAS_GCR && GCR_MODULE !== nothing
            # Use GCR.jl diagonal preconditioner
            try
                DiagonalPreconditioner = getfield(GCR_MODULE, :DiagonalPreconditioner)
                P = DiagonalPreconditioner(Float64)
                # Setup preconditioner
                setup_preconditioner! = getfield(GCR_MODULE, :setup_preconditioner!)
                setup_preconditioner!(P, matrix)
                return P
            catch e
                @warn "Failed to create GCR diagonal preconditioner: $e, using fallback"
                return Diagonal(1.0 ./ diag(matrix))
            end
        else
            # Simple diagonal scaling for fallback
            diag_vals = diag(matrix)
            # Avoid division by zero
            safe_diag = [abs(d) > 1e-14 ? 1.0/d : 1.0 for d in diag_vals]
            return Diagonal(safe_diag)
        end
        
    elseif precond_type == :ilu
        if HAS_GCR && GCR_MODULE !== nothing
            # Use GCR.jl ILU preconditioner if available
            try
                IncompleteLUPreconditioner = getfield(GCR_MODULE, :IncompleteLUPreconditioner)
                P = IncompleteLUPreconditioner(Float64, 1e-3, 1e-1)  # drop_tol, pivot_tol
                setup_preconditioner! = getfield(GCR_MODULE, :setup_preconditioner!)
                setup_preconditioner!(P, matrix)
                return P
            catch e
                @warn "ILU preconditioner failed: $e, falling back to diagonal"
                return setup_preconditioner(matrix, :diagonal, backend)
            end
        else
            # Use built-in incomplete LU if available
            try
                # Check if IterativeSolvers has ilu function
                if isdefined(IterativeSolvers, :ilu)
                    return ilu(matrix, τ=0.01)
                else
                    @warn "ILU not available in IterativeSolvers, using diagonal"
                    return setup_preconditioner(matrix, :diagonal, backend)
                end
            catch e
                @warn "ILU failed: $e, using diagonal preconditioner"
                return setup_preconditioner(matrix, :diagonal, backend)
            end
        end
        
    elseif precond_type == :amg
        # Algebraic multigrid preconditioner (future enhancement)
        @warn "AMG preconditioner not yet implemented, using diagonal"
        return setup_preconditioner(matrix, :diagonal, backend)
        
    else
        @warn "Unknown preconditioner type: $precond_type, using diagonal"
        return setup_preconditioner(matrix, :diagonal, backend)
    end
end

"""
    solve_poisson_fallback!(solver::Union{GCRPoissonSolver, FallbackPoissonSolver}, 
                           rhs_interior::Vector{Float64}, 
                           solution_field::Array{Float64}) -> NamedTuple

Emergency fallback solver when primary methods fail.
"""
function solve_poisson_fallback!(solver::Union{GCRPoissonSolver, FallbackPoissonSolver}, 
                                 rhs_interior::Vector{Float64}, 
                                 solution_field::Array{Float64})
    t_start = time()
    
    try
        # Simple conjugate gradient with minimal setup
        x0 = zeros(length(rhs_interior))
        A = solver.matrix
        b = rhs_interior
        
        # Manual CG implementation for maximum compatibility
        r = b - A * x0
        p = copy(r)
        x = copy(x0)
        
        rsold = dot(r, r)
        residual_norms = [sqrt(rsold)]
        
        maxiter = min(500, 2 * length(b))
        tol = 1e-6
        
        for iter in 1:maxiter
            Ap = A * p
            α = rsold / dot(p, Ap)
            x .+= α .* p
            r .-= α .* Ap
            rsnew = dot(r, r)
            
            residual_norm = sqrt(rsnew)
            push!(residual_norms, residual_norm)
            
            if residual_norm < tol
                scatter_interior_values!(solution_field, x, solver.interior_indices)
                solve_time = time() - t_start
                
                return (
                    converged = true,
                    iterations = iter,
                    residual_norm = residual_norm,
                    solve_time = solve_time,
                    solver_info = "Emergency CG"
                )
            end
            
            β = rsnew / rsold  
            p .= r .+ β .* p
            rsold = rsnew
        end
        
        # Did not converge
        scatter_interior_values!(solution_field, x, solver.interior_indices)
        solve_time = time() - t_start
        
        return (
            converged = false,
            iterations = maxiter,
            residual_norm = sqrt(rsold),
            solve_time = solve_time,
            solver_info = "Emergency CG (not converged)"
        )
        
    catch e
        @error "Even fallback solver failed: $e"
        solve_time = time() - t_start
        
        return (
            converged = false,
            iterations = 0,
            residual_norm = Inf,
            solve_time = solve_time,
            solver_info = "Total failure"
        )
    end
end

"""
    benchmark_poisson_solvers(domain::LShapeDomain{D}, h::Float64; 
                             backends=[:julia], rhs_type=:random) where D

Benchmark different Poisson solver configurations.

# Arguments
- `domain::LShapeDomain{D}`: Test domain
- `h::Float64`: Grid spacing
- `backends::Vector{Symbol}`: Backends to test
- `rhs_type::Symbol`: Type of right-hand side (:random, :manufactured)

# Returns  
- Performance comparison results
"""
function benchmark_poisson_solvers(domain::LShapeDomain{D}, h::Float64; 
                                  backends=[:julia], rhs_type=:random) where D
    
    println("Benchmarking Poisson Solvers ($(D)D)")
    println("="^50)
    
    # Generate test problem
    if D == 2
        nx, ny = size(domain.fluid_mask)
        rhs_field = zeros(nx, ny)
        solution_field = zeros(nx, ny)
        
        if rhs_type == :random
            rhs_field[domain.interior_mask] = randn(sum(domain.interior_mask))
        elseif rhs_type == :manufactured
            # Use manufactured solution for accuracy test
            for j in 1:ny, i in 1:nx
                if domain.interior_mask[i, j]
                    x_coord = -1.0 + (i-1) * h
                    y_coord = -1.0 + (j-1) * h
                    rhs_field[i, j] = -2π^2 * sin(π*x_coord) * sin(π*y_coord)
                end
            end
        end
    else  # 3D
        nx, ny, nz = size(domain.fluid_mask)
        rhs_field = zeros(nx, ny, nz)  
        solution_field = zeros(nx, ny, nz)
        
        if rhs_type == :random
            rhs_field[domain.interior_mask] = randn(sum(domain.interior_mask))
        end
    end
    
    results = []
    
    for backend in backends
        try
            options = LidDrivenOptions(solver=backend, verbose=false)
            solver = create_poisson_solver(domain, h, options)
            
            # Warmup run
            solve_poisson!(solver, rhs_field, copy(solution_field))
            
            # Benchmark run
            sol_test = copy(solution_field)
            result = solve_poisson!(solver, rhs_field, sol_test)
            
            push!(results, (backend=backend, result=result, solution=sol_test))
            
            println(@sprintf("%-15s: %s in %d iters, %.3f sec", 
                           string(backend),
                           result.converged ? "✓" : "✗",
                           result.iterations,
                           result.solve_time))
                           
        catch e
            println(@sprintf("%-15s: FAILED - %s", string(backend), string(e)))
        end
    end
    
    return results
end

"""
    determine_preconditioner_candidates(matrix::SparseMatrixCSC{Float64, Int}, backend::Symbol) -> Vector{Symbol}

Determine suitable preconditioner candidates based on matrix properties and backend.
"""
function determine_preconditioner_candidates(matrix::SparseMatrixCSC{Float64, Int}, backend::Symbol)
    n = size(matrix, 1)
    nnz_density = nnz(matrix) / n^2
    
    # Base candidates always available
    candidates = [:diagonal]
    
    # Add more sophisticated preconditioners based on problem size and density
    if n < 10000  # Small to medium problems
        push!(candidates, :ilu)
        
        if backend == :petsc
            push!(candidates, :gamg)  # Algebraic multigrid via PETSc
        end
    elseif n < 50000  # Medium problems
        push!(candidates, :ilu)
        
        if backend == :petsc || (HAS_GCR && GCR_MODULE !== nothing)
            push!(candidates, :gamg)
        end
    else  # Large problems
        # For very large problems, avoid expensive preconditioners
        if backend == :petsc
            push!(candidates, :gamg)
            push!(candidates, :hypre)  # HYPRE BoomerAMG
        end
    end
    
    # Check matrix properties for additional recommendations
    if is_symmetric_check(matrix)
        # For symmetric matrices, we could add specialized preconditioners
        if n < 5000 && backend != :julia
            push!(candidates, :cholesky)
        end
    end
    
    return unique(candidates)
end

"""
    is_symmetric_check(matrix::SparseMatrixCSC{Float64, Int}) -> Bool

Quick symmetry check for matrix (sample-based for large matrices).
"""
function is_symmetric_check(matrix::SparseMatrixCSC{Float64, Int})
    n = size(matrix, 1)
    if n <= 1000
        return issymmetric(matrix)
    else
        # Sample-based check for large matrices
        sample_size = min(100, n)
        indices = randperm(n)[1:sample_size]
        
        for i in indices
            for j in indices
                if abs(matrix[i,j] - matrix[j,i]) > 1e-12
                    return false
                end
            end
        end
        return true
    end
end

"""
    solve_poisson!(solver::AdaptivePoissonSolver, rhs_field::Array{Float64}, 
                  solution_field::Array{Float64}) -> NamedTuple

Solve Poisson equation with adaptive preconditioner selection.
"""
function solve_poisson!(solver::AdaptivePoissonSolver, rhs_field::Array{Float64}, 
                       solution_field::Array{Float64})
    
    best_result = nothing
    best_performance = Inf
    
    # If we have a current best, try it first
    if solver.current_solver !== nothing
        try
            result = solve_poisson!(solver.current_solver, rhs_field, copy(solution_field))
            
            if result.converged
                performance_score = compute_performance_score(result)
                
                # Update performance history
                push!(solver.performance_history[solver.best_preconditioner], performance_score)
                
                # Apply iterative refinement if enabled and needed
                if solver.refinement_enabled && result.residual_norm > solver.solver_options.tol * 0.1
                    refined_result = apply_iterative_refinement!(solver, rhs_field, solution_field, result)
                    if refined_result !== nothing
                        result = refined_result
                    end
                end
                
                scatter_interior_values!(solution_field, 
                                       extract_interior_values(solution_field, solver.interior_indices),
                                       solver.interior_indices)
                
                return result
            end
        catch e
            @debug "Current best solver failed: $e"
        end
    end
    
    # Try all candidates and find the best one
    for (i, precond) in enumerate(solver.preconditioner_candidates)
        try
            # Create temporary solver with this preconditioner
            temp_solver = GCRPoissonSolver(solver.matrix, solver.id_map, solver.interior_indices,
                                         backend = solver.backend,
                                         preconditioner = precond,
                                         tol = solver.solver_options.tol,
                                         maxiter = solver.solver_options.maxiter,
                                         restart = solver.solver_options.restart,
                                         verbose = false)
            
            test_solution = copy(solution_field)
            result = solve_poisson!(temp_solver, rhs_field, test_solution)
            
            if result.converged
                performance_score = compute_performance_score(result)
                
                # Update performance history
                push!(solver.performance_history[precond], performance_score)
                
                if performance_score < best_performance
                    best_performance = performance_score
                    best_result = result
                    solver.best_preconditioner = precond
                    solver.current_solver = temp_solver
                    
                    # Copy the best solution
                    solution_field .= test_solution
                end
                
                # Early termination if we find a very good solution
                if result.residual_norm < solver.solver_options.tol * 0.01 && result.iterations < 20
                    break
                end
            end
            
        catch e
            @debug "Preconditioner $precond failed: $e"
            push!(solver.performance_history[precond], Inf)
        end
    end
    
    if best_result !== nothing
        # Apply iterative refinement if enabled
        if solver.refinement_enabled && best_result.residual_norm > solver.solver_options.tol * 0.1
            refined_result = apply_iterative_refinement!(solver, rhs_field, solution_field, best_result)
            if refined_result !== nothing
                best_result = refined_result
            end
        end
        
        return best_result
    else
        # All methods failed, fall back to emergency solver
        @warn "All adaptive preconditioners failed, using emergency fallback"
        return solve_poisson_fallback!(solver, 
                                     extract_interior_values(rhs_field, solver.interior_indices), 
                                     solution_field)
    end
end

"""
    compute_performance_score(result) -> Float64

Compute a performance score combining solve time, iterations, and convergence quality.
Lower scores are better.
"""
function compute_performance_score(result)
    if !result.converged
        return Inf
    end
    
    # Weighted combination of time, iterations, and residual quality
    time_weight = 1.0
    iter_weight = 0.01  # Each iteration adds 0.01 to score
    residual_weight = 100.0  # Penalty for poor convergence
    
    score = time_weight * result.solve_time + 
            iter_weight * result.iterations + 
            residual_weight * max(0, log10(result.residual_norm + 1e-16))
    
    return score
end

"""
    apply_iterative_refinement!(solver::AdaptivePoissonSolver, rhs_field::Array{Float64},
                               solution_field::Array{Float64}, initial_result) -> Union{NamedTuple, Nothing}

Apply iterative refinement to improve solution accuracy.
"""
function apply_iterative_refinement!(solver::AdaptivePoissonSolver, rhs_field::Array{Float64},
                                    solution_field::Array{Float64}, initial_result)
    
    if solver.current_solver === nothing
        return nothing
    end
    
    max_refinement_steps = 3
    refinement_tol = solver.solver_options.tol * 0.1
    
    current_solution = copy(solution_field)
    rhs_interior = extract_interior_values(rhs_field, solver.interior_indices)
    
    for step in 1:max_refinement_steps
        # Compute residual: r = b - A*x
        x_interior = extract_interior_values(current_solution, solver.interior_indices)
        residual = rhs_interior - solver.matrix * x_interior
        residual_norm = norm(residual)
        
        if residual_norm < refinement_tol
            # Refinement successful
            refined_result = (
                converged = true,
                iterations = initial_result.iterations + step,
                residual_norm = residual_norm,
                solve_time = initial_result.solve_time,  # Approximate
                solver_info = "$(initial_result.solver_info) + refinement($step)"
            )
            
            solution_field .= current_solution
            return refined_result
        end
        
        # Solve correction: A * delta_x = residual
        correction_field = zeros(size(solution_field))
        scatter_interior_values!(correction_field, residual, solver.interior_indices)
        
        delta_solution = zeros(size(solution_field))
        
        try
            correction_result = solve_poisson!(solver.current_solver, correction_field, delta_solution)
            
            if correction_result.converged
                # Apply correction with damping
                damping = 0.8  # Damping factor to ensure stability
                current_solution .+= damping .* delta_solution
            else
                break  # Refinement failed
            end
        catch
            break  # Refinement failed
        end
    end
    
    return nothing  # Refinement did not improve the solution significantly
end

"""
    get_adaptive_solver_stats(solver::AdaptivePoissonSolver) -> Dict

Get performance statistics for the adaptive solver.
"""
function get_adaptive_solver_stats(solver::AdaptivePoissonSolver)
    stats = Dict{String, Any}()
    
    stats["best_preconditioner"] = solver.best_preconditioner
    stats["candidates_tested"] = collect(keys(solver.performance_history))
    
    # Compute average performance for each preconditioner
    avg_performance = Dict{Symbol, Float64}()
    for (precond, times) in solver.performance_history
        if !isempty(times) && !all(isinf, times)
            valid_times = filter(t -> !isinf(t), times)
            avg_performance[precond] = isempty(valid_times) ? Inf : mean(valid_times)
        else
            avg_performance[precond] = Inf
        end
    end
    
    stats["average_performance"] = avg_performance
    stats["refinement_enabled"] = solver.refinement_enabled
    stats["total_solves"] = sum(length(times) for times in values(solver.performance_history))
    
    return stats
end